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Abstract 

Extreme value theory for chaotic deterministic dynamical systems is a rapidly expanding 

00 . 

CN| , area of research. Given a system and a real function (observable) defined on its phase space, 

extreme value theory studies the limit probabilistic laws obeyed by large values attained by 

L—\ • the observable along orbits of the system. Based on this theory, the so-called block maximum 

method is often used in applications for statistical prediction of large value occurrences. In 

a . 

this method, one performs statistical inference for the parameters of the Generalised Extreme 
Value (GEV) distribution, using maxima over blocks of regularly sampled observations along 
an orbit of the system. The observables studied so far in the theory are expressed as functions 
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\^ , of the distance with respect to a point, which is assumed to be a density point of the system's 
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invariant measure. However, this is not the structure of the observables typically encountered 
in physical applications, such as windspeed or vorticity in atmospheric models. In this paper we 
consider extreme value limit laws for observables which are not functions of the distance from 
a density point of the dynamical system. In such cases, the limit laws are no longer determined 
by the functional form of the observable and the dimension of the invariant measure: they also 
depend on the specific geometry of the underlying attractor and of the observable's level sets. 
We present a collection of analytical and numerical results, starting with a toral hyperbolic 
automorphism as a simple template to illustrate the main ideas. We then formulate our main 
results for a uniformly hyperbolic system, the solenoid map. We also discuss non-uniformly 
hyperbolic examples of maps (Henon and Lozi maps) and of flows (the Lorenz63 and Lorenz84 
models). Our purpose is to outline the main ideas and to highlight several serious problems 
found in the numerical estimation of the limit laws. 
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I. INTRODUCTION 



Background on extreme value theory 



Classic extreme value theory concerns the probability distribution of unlikely (large) 



events, see 



aa 



10 



13 



24 



37 



43| . Given a stochastic process Xi,X 2 ,... governed 
by independent identically distributed random variables, let M n be the random variable 
defined as the maximum over the first n occurrences: 

M n = max(Xi,. . . ,X n ). 

This variable has a degenerate limit as n — > oo, and therefore it is necessary to consider 
a rescaling. Suppose that there exist sequences a n > and b n EM. such that the rescaled 
variable a n (M n — b n ) converges to a non-degenerate distribution. That is 

lim P (a n {M n - b n ) <x) = G(x) (1) 

Then extreme values theory asserts that the limit G(x) can only be one of three different 
types: the Gumbel, Weibull and Frechet parametric families of probability distributions. 
These three families can be combined into a single three-parameter family having distri- 
bution function 

x — n 



G(x) = exp 



£ 



i -ye 



(2) 



a 

defined on the set |x| 1 + £ (— ^-^-) > 0}, where the parameters satisfy — oo < [i < oo, 
a > and — oo < £ < oo. Eq. ([2]) is called the generalised extreme value (GEV) family of 
distributions. The subset of the GEV family with £ = is interpreted as the limit of ([2]) 
as £ — y 0, leading to the Gumbel family (with parameters /i and a). 

In the applications the GEV family is particularly useful to predict the probability of 
occurrence of future large values of a quantity, given a sample of past experimental mea- 
surements of that quantity. The so-called block maximum method is frequently used in 
this setting. Here one extracts a sub-sample of maxima over data blocks: in environ- 
mental and climate contexts one often uses blocks of length one year, hence the name of 
annual maximum method. One then estimates the parameters (/i, a, £), assuming that 
that the block maxima form a random sample drawn from a GEV distribution with un- 



known parameters. Maximum likelihood is a common estimation method [10|: in this 
case, standard asymptotic theory also provides confidence intervals (uncertainties) for the 
point estimates. The estimated GEV parameters and associated uncertainties can then 



be used to derive other quantities of interest, such as return periods for given return levels 



of the variable of interest, see the above references and 



16 



17 



52 



53| for examples. 



Extremes in deterministic systems 



Recent work has extended the domain of extreme value theory to the setting of chaotic 



deterministic dynamical systems 



BQ 



18 



22 



23, m, |3l|. We briefly outline the dif- 



ference of our problem setting as opposed to the above results. Suppose that we have a 
dynamical system (X, u, f), where X is a d- dimensional Riemannian manifold, / : X — > X 
a measurable map and v an /-invariant probability measure. Assume that there is a com- 
pact invariant set A C X which supports the measure v. Specifically, our main interest 
is the situation where A is a strange attractor and v is a Sinai-Ruelle-Bowen (SRB) mea- 
sure [59]. Given an observable <p : X — > M U {+00} we study extreme value limit laws for 
the stationary stochastic process Xi, X 2 , ■ ■ ■ defined by 



Xi 



of 



i-1 



I > 1. 



(3) 



The theoretical work cited above focused on the special case where </> has the form 



4>{p) = 0(dist(p,pAf)), P6^, 



(4) 



where g : [0, +00) — y R is a measurable function of the distance dist(-, •) in X and pu is a 
density point of v. However, typical observation functions used i n ap plications are not of 
this form. Consider for example the quasi-geostrophic model of 16|, ll7|]: this model was 



conjectured in [40| to possess a compact (bounded) strange attractor in its (unbounded) 



] Q 



56[ are the system's total energy, the 



phase space. The observables used in [16|, H7J, [52J, 

wind speed and vorticity at a gridpoint in the lower level. These observables can be 

written as 



<Mp) =p T Ep, <fiw(p) = \\Wp\\, (j> v {p) = Vp, respectively, 



(5) 



where p is a point in the phase space X = M. d , \\ ■ \\ denotes the Euclidean norm and 
E G R dxd , W G R 2xd , V G R lxd are matrices. None of the observables in (jSJ) has the 
form (HI). In fact this situation is to be expected in many, if not in most, observables 

nri 

found in applications, including the atmospheric and oceanic models of [7J, |47J. Although 
observables such as flH]) are usually unbounded in the system's phase space, the system's 
attractor A is usually bounded due to the presence of dissipative processes in the models. 
Therefore, time series of such observables should be expected to have an upper bound 



and, hence, large values typically obey Weibull limit distributions, see 16|, |52| for a more 
detailed discussion. 



Sketch of the results 

In this paper we will consider observables which are not a function of the distance from 
a point pm as in (j4]). Such observables include cases, like those just mentioned, where 
has no upper bound in the phase space, although time series of on the system's attractor 
are bounded. Hence we will restrict to the Weibull case in our numerical examples. For 
comparison with the already available theory, we also consider cases when is maximised 
at a point pm, where however pu may or may not be a density point for the invariant 
measure v. 

In such situations, to determine the form of the limiting GEV distribution G{x) becomes a 
much more delicate problem: G(x) is no longer determined by the functional form of 0(p) 
and by the dimension of the SRB measure z/, but also critically depends on the geometry 
of the attractor A C X and of the level sets of 0. A careful analysis is required even 
if we assume that takes the form of (jlj), but allowing instead that pm ^ A. Without 
attempting an exhaustive analysis of all possible cases, we focus on selected examples to 
illustrate the key ideas of our approach, in view of applications to a given system and 
observable. 

To whet the Reader's appetite, we here anticipate one of the results of this paper. For 
/ we consider the solenoid map [6( embedded in M 3 = {(x, y,z)}\ this system possesses 



a strange attractor A which is locally the product of a Cantor set with an interval 301 ]. 
where the interval represent the leaves of a one-dimensional unstable manifold W u . For 
the observable we take 0(x, y, z) = ax + by + c + d, which is clearly not of the form (J4]): 
rather, resembles the vorticity observable defined in ([5]). For this pair of system and 
observable we obtain the formula 

-\ = Y + d ' (6) 

for the tail index £ of the limiting GEV distribution. Here d u = dim(W u ) = 1 and d s is 
the dimension in the stable direction, which in this case is given by dim^(A) — 1, where 
dim# denotes the Hausdorff dimension. Loosely speaking, the factor 1/2 in ([6]) is obtained 
under a generic condition of quadratic tangency between a local unstable manifold within 
the attractor A and the level sets of the observable 0. As we shall argue, we believe 



formula (]6]) to be valid, or at least sufficiently informative, for a large class of pairs (/, 0) 
of systems-observables. However, we also discuss examples where formula ([6]) has to 
be modified to take into account the local geometry or the local scaling of the invariant 
measure of the attractor, or the local behaviour of the level sets of the observable. We will 
restrict our discussion to the tail index £, which is usually the most delicate parameter to 
estimate in the analysis of extreme values: see e.g. |2|, [l0[ and 15( for the link between the 



normalising constants a n , b n and the other two GEV parameters ft, a . We note, however, 
that our numerical procedure also provides estimates for the latter two parameters, see 
Appendix [A] 

Outline of the paper 

This paper is organised as follows. The general framework of extremes in dynamical 
systems is presented in Section [Til Our main theoretical results are formulated in Sec- 
tions IIHI to [VI] for specific dynamical systems. As a particularly simple example, we 
consider Thorn's map in Section IIHI to illustrate the main ideas in our approach. Then 
in Section [TV] we discuss the solenoid attractor, which displays several features which 
are found in many concrete physical systems. Section |V] presents results for two non- 
uniformly hyperbolic systems, the Henon and Lozi maps. In Section [VI] we examine two 



prototypical flows with chaotic dynamics due to Lorenz [38|, [39|. In all of the sections, 
analytical calculations precede numerical simulations, where the latter aim to show the 
typical behaviour and estimation problems which can be expected to occur. We return in 
Section IVIII to the general relevance of our approach. 

II. EXTREMES IN DYNAMICAL SYSTEMS: THE GENERAL PROBLEM 
SETTING 

We consider a measure preserving system (X, u, /) with a compact attracting set A C X. 
Given an observable : X — y M and a threshold u G R, we define the level regions L + (u) 
(resp. level sets L{u)) as follows: 

I+{u) = {p e X : <f>{p) > u}, L(u) = {peX :(f)(p)=u}. (7) 

For the reasons discussed in the Introduction, we consider observables which achieve a 
finite maximum within A, although the observable themselves could be unbounded in X. 



We define 



u = sup0(jo). 



Since A is compact there exists (at least) one point p£ A for which cf)(p) = u. We will 
assume that such an extremal point p is unique. Given our focus on the Weibull case 
(again, see the Introduction) we consider sequences u n = u/a n + b n for which the limit 

lim nv(L + (u n )) := r(«) (9) 

n 

exists. From the theory of [36|, we can choose b n = u, and we take a n — > oo. The precise 
form of a n depends on the regularity of 0, and the regularity of the density of v in the 
vicinity of the extremal point p. In general a n will be a power law in n, and r(u) will 
be regularly varying in u. If we now consider the process M n = max(Xi, . . . ,X n ) with 
X n — (j) o f n , then we investigate to what extent the following statement is true: 

nv{p : 0(p) > u n } -> r{u) <£> v{M n < u n } -> e~ r{u \ (10) 

If t(u) = u a , then the process M n is described by a GEV distribution with tail index 
£ = — 1/a. The statement ( TTOl) is shown to hold for a wide class of dynamical systems, such 
as those governed by non-uniformly expanding maps, and systems with (non)-uniformly 
hyperbolic attractors, [9j, |28[ . However the current theory for analysing extremes in dy- 
namical systems assumes that the level regions L + (u) introduced in ([7j) are described by 
balls, and moreover that these balls are centred on points in A that are generic for v. 
These assumptions allow for the tail index to be expressed in terms of local dimension 
formulae for measures. 

In this article we do not assume that the level sets are balls: for example we consider 
observables of the form <j)(p) = 4>(xi, . . . ,Xd) = Y2i \ x i\ ai i where the level sets have cusps 
or are non-conformal. We also consider observables 4>(p) = ^^CjXj, for which the level 
sets are hyperplanes. For observables of these types (also compare with (jSJ)) the standard 
machinery does not immediately apply. The first problem is to determine the sequence 
u n and the limit t(u) defined in Q. Even if the measure v is sufficiently regular then 
the sequence u n will depend on the geometry of the attractor close to where <p(p) achieves 
its maximum value on A, in addition to depending on the form of 0. In Section IHII 
we illustrate the various geometrical scenarios that can arise using an hyperbolic toral 
automorphism as a simple example. When v is a more general SRB measure, then even for 
uniformly hyperbolic systems (such as the solenoid map) it becomes a non-trivial problem 
to determine u n and r(u). We discuss this scenario in Section HVl The second problem is 
to verify statement ( fTUl) for the class of observables under consideration. This relies upon 



checking two conditions D 2 (u n ) , D' (u n ) , see [19|, |28[. We summarise these conditions as 



follows. For integers £, / let M tj i = max{X t+ i, X t +2, . . . , X t+ i}, with M j := Mi. Then: 

(D2(u n )) We say condition D2{u n ) holds for the process X , Xi, . . . , if for any integers l,t 
and n we have 

\v(Xi > u n ,M t) i < u n ) - u(X 1 > Un)v(Mi <Un)\ < j(n,t), 

where 7(71, t) is non-increasing in t for each n and 727(72, t n ) — > as n — > 00 for some 
sequence t n = o(n), t n —¥ 00. 

(D'(w n )) We say condition D (%) holds for the process Xi, X2, ■ ■ ■ , if 

[n/fc] 

lim limsupn N, ^(^1 > u n,Xj > u n ) = 0. (11) 

If the level sets have complicated geometry, or if the measure v is supported on a fractal 

ic systems, 



28|. In this 



set then these conditions must be carefully checked. For uniformly hyperbo 
and for observations that are functions of balls these conditions are checked in 
article we consider the computation of the GEV tail index £ for more general observations 
and contrast to results known for observations that are functions of balls. We focus on 
particular examples to highlight how the geometrical features of the level sets and the 
attractor feed into the computation of the tail index £, without attempting an exhaustive 
analysis of all possible cases. We also discuss the computation of the tail index for non- 
uniformly hyperbolic systems such as the Lozi map and Henon map, and also for Lorenz 
flows, again for general observations. 

III. A PROTOTYPICAL EXAMPLE: THOM'S MAP 

Let T 2 = IR 2 mod 1 be the 2 dimensional torus. Thorn's map / : T 2 — > T 2 , also known 
as Arnold's cat map, is the hyperbolic toral automorphism defined by 

f(x,y) = (2x + y,x + y) modi. (12) 

This system is Anosov and it has Lebesgue measure v on the torus T 2 as the (unique) 
invariant measure. With this example we want to study the role of the observable in 
determining extreme value laws. For this purpose we will consider / as a map of IR 2 
having the square [0, l) 2 as the invariant set. In other words, X = IR 2 and A = [0, l) 2 , 
hence A is not an attractor, strictly speaking. The advantage is that this allows us to take 



functions of M 2 as observables, rather than functions of T 2 . In this way, we can construct 
observables which are maximised at points in the interior or in the complement of A and 
whose level sets have different shapes. 

The main point of this section is that the value of the tail index is determined by the 
interaction between the shape of level sets (E]) of the observable and the shape of the 
support of the invariant measure (colloquially, the geometry of the attractor). To illustrate 
our ideas, and without attempting to cover all possible cases, we consider the following 
two observables <f) a , <p a ^ : M 2 — > M. 

<p a (x,y) = 1 - distO,p M ) a , with p = (x,y) E R 2 . (13) 

<f>a,b(%, y) = l-\x- x M \ a -\x- y M \ b , (14) 

where, given our focus on the Weibull case, we require a,b,a > 0. Both observables are 
maximised at a point pu = (xm,Vm) £ R 2 - When pu is in the interior of A, observ- 
able (TTBl) has the form so far analysed in the mathematical literature about extremes in 
dynamical systems, but we will also consider the case pm $■ A. Observable (JHj) has been 
chosen to illustrate the effect of the shape of the level sets: the level regions of ([14]) are 
not balls unless a = b = 2, in which case (TT4|) can be written as (Tl3|) for a = 2. The 
next subsection contains our analytical results, numerical simulations are postponed to 
Section ME 

A. Analytical calculations 

The level regions L + (u) as defined in ([7]) are always balls for observable ffTBl . However 
three (main) different situations occur, depending on the location of the point pm relative 
to the support of the invariant measure, see the sketch in Figure [TJ 

Theorem III.l. Let £ be the tail index of the GEV limit distribution associated to the 
process M n = max(Xi, . . . , X n ) with X n = <po f n ! where f is the map [W\l and <p a : M 2 — > 
IR is the observable in (T73p- Then for v-a.e. pu = (xm,Vm) E K 2 we have: 

forp M eA; (15) 



1 2 

I a 
1 3 

1 

-r 2 



for p M g" A, with either y AI G (0, 1) or xm £ (0, 1); (16) 

for p M g A, with both x M , yu & [0, 1]; (17) 



For observable ( !T4l) the shape of the level sets L(u) depends on a and b. For example, L(u) 
has a convex elliptic-like shape when both a, b > 1 (see the sketch in Figure [2] (A)), or 



an asteroid-like shape when both a,b < 1, (see Figure |2](B)). Clearly various possibilities 
arise, depending on the geometry of the level sets, on whether the point pm is in the 
interior of A and on the local geometry of A near the extremal point p = (x, y) with 
minimum distance from pm- 

Theorem III. 2. Let £ be the tail index of the GEV limit distribution associated to the 
process M n = max(Xi, . . . ,X n ) with X n = <fi o f n ! where f is the map (TB) and <p a ^ '■ 
R 2 — > K is the observable in [Lj\) . Then for v-a.e. pu = (xm,Vm) G K 2 we have: 

-t = - + t forp M eA. (18) 

4 a b 

To prove Theorems IIII.lHIH."2l the main step is to determine the explicit sequence u n and 
functional form of t(u) as defined in (El). We will not give the verification of D2(u n ), 



D'(u n ) as this follows step by step from 28j for this class of observables. The main proof 
of Theorem IIII. II is contained in Lemmas |III.3|III.4I and IIII.51 The proof of Theorem IIII. 21 
is given in Lemma 1111.61 We do not further discuss case (C) of Figure |2j or the other 
configurations not covered by Figure [2j 

Lemma III. 3. Suppose pu e int(A) = (0, l) 2 and <fi takes the form of (T73)|, then £ = 

-a/2. 

Proof. If pm £ int(A) (see Figured] (A)), then we see that 

nv{<f)(x, y) > u n } = nv{x < (1 - u n ) 1/a } 

(19) 

= n{l-u n f/ a . 

Thus the correct scaling laws are a n = n a < 2 , b n = 1 and r(u) = {—u) 2 ' a . D 

Lemma III. 4. Suppose pu ^ A = [0, l] 2 and <fi takes the form of |73|). If Vm G (0, 1) or 
x M e (0, 1) then £ = -2/3. 

Proof. If pm $■ A then there will exist a unique extremal point p = (x, y) G A where 4>(p) 
achieves its supremum with value u as in (JSJ). Since yM G (0, 1) or xm G (0, 1) then this 
point p will not be a vertex of <9A, see Figured] (B). The scaling u n will be chosen to so 
that 

nu{(p(x,y) > u n } = nv{p = (x, y) G A : d(p,p M ) < (1 - u n ) 1/a } -> r{u). (20) 

The middle term is no longer n{\ — u n ) l l a since the level region that intersects A is not 
a ball. We first of all set u n = u/a n + u so that 

v{<i>{x,v) >u n } = ulp=(x,y) GA: (1 - u) l/a < d(p,p M ) < (l-u-^-j 

'(21) 



To choose a n we first note that the level set Liu 1 ' 01 ) as defined in (|7]) is a circle that 
is tangent to dA (since the extremal point p is not a vertex). However the level set 
L((u — ^-) 1 / a )) crosses dA transversally (and is concentric to L{u l ^ a )). Thus a basic 
geometrical calculation gives: 

v{cf>(x, y) > u n } = {(1 - u - u/a n fl a - (1 - u) 1 '^' 2 

l/a ^ 3 / 2 

a„(l — u) 



'l-uf/ 2a l (l J- r ) -1 



(22) 



(l-n) 3 / 2 * ^^T + O 



1 3/2 
U „ / 1 N 



cra„(l — u) V O; 



,2 



Setting a n = n 2//3 implies that r(n) = (— u) 3//2 . D 

Lemma III. 5. Suppose puj ^ A and <p takes the form of |73j). If both %, ?/m ^ [0, 1] 
then £ = -1/2. 

Proof. Without loss of generality we consider p M = (xm,Vm) i n the upper right hand 
quadrant as in Figured] (C). Also, we assume that xm = 1 + Acos6>, y M = 1 + A sin 6* for 
A > and 9 e (0,7r/2). For such values of (xm,Vm), the corner point (1,1) G dA will 
always maximise <fi. The proof is identical to Lemma IIII.4I except that the level sets are 
not tangent to dA at (1, 1), as illustrated in Figured] (C). A simple geometric argument 
shows that for a disk D(e) of radius A + e, centred at (xm, Vm) we have 

Leb{p = (x, y) e M n D(e)} = C(e 2 ), 

where Leb denotes the Lebesgue measure. Setting u n = u/a n + u and comparing to (|22|) . 
we obtain: 

2 



u{cj>(x, y)>u n } = {(l-u- u/a n y/ a - (1 - n) 1 ^}' 



l- n u ^ y /Q -i 

a n (l -u)J 
aa n {l-u) V a n 



(23) 



Setting a n = n 1 / 2 implies that t{u) = (—u) 2 . □ 

This concludes the proof of Theorem IIII.ll For the proof of Theorem IIII.2I we have the 
following lemma. 

Lemma III. 6. Suppose that pm £ int(A) = (0, l) 2 and takes the form of (T$ ). Then 
for u < 1 we have that Leb(L(u)) = C(l — u)* + b for some Cq < C < 4 where C$ > 0. 



Proof. Let u — 1 — £. For e sufficiently small the level region can be written as 

L + {u) = {(x, y) G int(A) : \x\ a + \y\ b < e}. (24) 

The area of this set is bounded from above by the area of a rectangle of sides 2e 1 ' a and 
2e x l h . Also, for any q G (0, 1), the area of the set is bounded from below by that of a 
rectangle of sides 2q l / a e l l a and 2(1 — q) l / b e l / b , so we can choose Cq = 4g 1 / a (l — q) l / h . □ 

Hence if (xm,Vm) G int(A), we see that 

nu{p = (x,y) : (f>(p) > u n } — > (— u)^ + b with a n = n^+z, b n = 1. (25) 

B. Numerical results 



As formula ( 1181) shows, one of the main ingredients in determining the tail index is the 
shape of the level sets of the observable. We here fix a = 2 and consider two values of b, 

is less 



10|, we 



namely 6 = 1 and b = 3.5. In both cases, the value of £ expected according to ( TT8I) 
than —0.5. Since the maximum likelihood estimator is not regular for £ < —0.5 
resort to the method of L-moments for the numerical estimation of the GEV parameters. 
See Appendix [A] for details on our procedure for parameter estimation and associated 
uncertainties. 

In Figure [3] we examine the sensitivity of the numerical estimates of £ with respect to the 
block length used to compute the maxima. Essentially no significant variations are found 
for block lengths larger than 1000. Hence, we fix Nuockien — 10000 and conduct a study 
of the dependence of the tail index on the parameter b of the observable ( fl8|) . Figure |4] 
shows a good agreement with the theoretical predictions of ffT8|) for a range of values of 
b. In this example the level regions of the observable of the form L + (u) as in ([7]) are fully 
contained in the interior of A = [0, l) 2 , at least for sufficiently large values of the threshold 
u. The only peculiarity is the non-circular shape of L + (u), see Figure |2] and compare with 
Lemma IIII.61 

We now consider a case where the level regions L + (u) are not fully contained in A. We 
take observable (Tl3|) and vary the location of the point pm = {xm,Vm)- Starting from 
values xm, Dm such that pu is in the interior of A we increase xm across 1, bringing p M 
in the region where yu G (0, 1) but xm $■ (0, 1). This transition is illustrated in panels 
(A) and (B) of Figured] and the two situations correspond to f fl5|) and ffT6l) respectively. 



Figure [5] shows the sensitivity of the numerical estimates of £ with respect to the block 
length used to compute the maxima for four values of xm ■ Convergence to the theoretical 
value ( fl5l) is achieved already with block lengths of a few hundred iterates when the point 
Pm is in the interior of A (panel (A), xm — 0.9) and on the boundary of A (panel (B), 
xm = 1-0). When pu is in the complement of A = [0, l] 2 but close to its boundary, 
then very large block lengths {Nj,i oc ki en > 10 5 ) are required to achieve convergence to 
the theoretical value of ( IT6|) (panel (C), xm = 1-01). When pu is further away from 
[0, l] 2 shorter block lengths of about 10 4 iterates already guarantee convergence to the 
theoretical value of ( TT6|) (panel (D), xm = 1-1). 

Figure [6] (A) shows the discontinuity of £ in the transition between the situations of panels 
(A) and (B) in Figure [Q The figure shows the estimated value of £ as a function of xm 
where x/m is kept constant and a fixed block length is used. As the point pu exits A, 
the value of £ has a jump from the value of (fT5j) to that of (fl6|) . However, the numerical 
estimation does not resolve this jump unless large block lengths (Nuockien > 10 5 iterates) 
are used. 

Lastly, Figure [6] shows the discontinuity of £ from (Tl6l) to ( fl7l) . at the transition between 
the situations of panels (B) and (C) in Figured! This transition is not resolved accurately 
even with block lengths of 10 5 . From the numerical point of view, very large block lengths 
are required near the transition to detect the change of scaling between ff22l) and fl23|) . 

The example discussed in this section is admittedly somewhat artificial. It has been chosen 
to clearly illustrate the main ideas and the problems which are found in the numerical 
estimation, without the additional complications due to higher dimensionality of phase 
space and fractal nature of the attractors. In the next section, we consider a situation 
which is closer to what one can expect in concrete physical systems. 

IV. UNIFORMLY HYPERBOLIC ATTRACTORS: THE SOLENOID MAP 

Consider the solid torus as the product of T = M/Z times the unit disc in the complex 
plane B# = {z G C| \z\ < 1}. Then the solenoid map is defined as follows: 

/ A :TxD4 TxD 

(26) 
(ip,w) h-> (2^, Xw + Ke l2 ^) . 

In order to have the map well defined we need K + XR < R and XR < K. For 
our purposes it is convenient to have the torus embedded in M 3 . Consider Cartesian 



coordinates (x, y, z) e R 3 and define corresponding cylindrical coordinates r, ip, z by 
x = rcos(^) and y = rsm(ip). Then the torus of width R can be identified with the 
set D — {(r — l) 2 + z 2 < -R 2 } for i? < 1. The torus T x © R (with coordinates (ip, u + iv)) 
can be identified with D taking r = 1 + u and 2; = v. We thus obtain an embedded 
solenoid map 

g x :D^D, g x (^,r, z) = (2^, 1 + iTcos(V') + X(r - 1), Ksin(^) + Az). (27) 

The solenoid attractor is defined as the attracting set of the map g x : 

For A < I we have 

dim^(A) = l + -^-, (28) 



ogA" 1 



45| . We consider the following observables 



where dim// denotes the Hausdorff dimension 
(j> a , <j) abcd : M 3 -> R: 

(j) a (x,y,z) = l-dist(p,p M ) a , with p = (x,y,z) eM 3 , (29) 

(pabcd(x, y, z) = ax + by + cz + d, (30) 

Observable (1291) is maximised at a point pjy/ G M 3 , whereas (130]) is unbounded in the 
phase space R 3 (except for the trivial choice a = b = c = 0). Note that the vorticity 
observable <pv of (|5]) has the same general form as ( 1301) . Our theoretical expectations are 
first discussed in Section IIV Al followed by numerical results in Section IIVBI 

A. Analytical calculations 

For observables which are functions of distance we have the following result. It is not 



explicitly stated in the literature but the proof follows straightforwardly from 



28 



Theorem IV. 1. Let £ be the tail index of the GEV limit distribution associated to the 



process M n = max(ATi, . . . , X n ) with X n = 4>og x , where g x is the map pT7| ) and (f) a : R — > 
R the observable of / T^Pj) . where pu £ A. Then we have: 

1 dim// (A) 



i a 



(31) 



More interesting considerations arise for the observable (J30l) . As a simple case, consider 
first the degenerate solenoid with A = and take a planar observable := ax + by + d, 



thus reducing the problem to the (x, y)-plane. In this case we have the trivial dimension 
formula dim// (A) = 1 since A is a circle. However, for computing the tail index we lose 
a factor of 1/2 due to the geometry of the level set. Indeed, level sets are straight lines 
within the (x, y)-plane, and at the extremal point p = (x,y) the critical level set L(u) is 
tangent to A. Since the tangency is quadratic, we find that 

u(L + (u - e)) = m 7U { 7 "(p) n L + {u - e)} = 0(y/e). (32) 

Here 7 u (p) is the unstable manifold through p (i.e. it is the unit circle), and m 7 « is the 
one-dimensional conditional (Lebesgue) measure on 7 u (p). Hence 

— = dim// (A) = -. 

f V ; 2 2 

The mechanism described above is similar to that illustrated for Thorn's map in Fig- 
ure H](B), leading to formula ( 1T6|) : indeed, there we have dim// (A) = 2, yielding the value 
3/2 for the tail -l/£. 

For A > 0, the attractor has more complicated geometry and is locally the product of a 



Cantor set with an interval 30]. Planar cross sections that intersect A transversely form a 
Cantor set of dimension dim// (A) — 1 = — log 2/ log A. To calculate u(L + (u — e)) we would 
like to repeat the calculation above using equation ( l32l) . but now the set of unstable leaves 
that intersect L + (u — e) form a Cantor set (for each e > 0). The extremal point p where 
<j)(p) attains its maximum on A forms a tip of A relative to L(u). Such a tip corresponds 
to a point on p G A whose unstable segment j u (p) is tangent to L(u) at p, and moreover 
normal to V0(p) at p. Given e > 0, we (typically) expect to find a Cantor set of values 
t G [0, e] for which the level sets L(u — t) are tangent to some unstable segment 7" C A. 
For other values of t these level sets cross the attractor transversally. Given (fixed) eo > 
we can define the tip set T = T(eo) C A as follows: let T P 7 u (p) be the tangent space to 
7" at p. Then we define 

r = {p G L + (u - e ) D A : T pl u (p) ■ V0(p) = 0}. (33) 

This tip set plays a role in proving the following result, which in turn provides us with 
information on the form of the tail index £. 

Proposition IV. 2. Suppose that g\ is the map (2l\ ) and = (p a bcd- Define r(e) = 



v{(j){p) > u — e}. //dini//(r) < 1, then modulo a zero measure set of values (a,b,c,d), 
7"(e) is regularly varying with index 1/2 + dim//(r) as e — > 0. 

We give a proof below. Based on this proposition we conjecture that 

-I = di m „(A)-I = I + J^_. (34) 



We outline the main technical challenges that would need to be overcome to prove this 
conjecture. Firstly, conditions D{u n ) and D'(u n ) should be checked. We believe that this 



should follow from 28j, however the proof would be non-standard due to the level set 



geometry. Secondly, we would claim that dim#(r) = dim#(A) — 1. The proof of this 



would utilise the techniques used in 30] to analyse the regularity of the holonomy map 



between stable disks. In particular, the Authors of [30( show that the holonomy map is 
Lipschitz on a set of full dimension. However, it does not automatically follow that the 
holonomy map between T and A fl D is Lipschitz (for a disk D transverse of A), but we 
believe that it is for general planar observations. 

Proof of Proposition \IV.Si For each e < eo, consider the set T(e) C T fl L + (u — e). Then 
for each p G T(e), there exists t < e such that 7 u (p) is tangent to L(u— t). If the observable 
4> takes the form of ( !30l) . then by the same calculation as ( l32l) we obtain 



m 7 »{f(p)nL + (M-e)} = 0(v^i). (35) 

Thus to compute u(L + (u — e)), we integrate ( )35l) over all relevant t < e using the measure 
u-p, which is the measure v conditioned on T. Provided dim^(r) < 1, the projection of 
T onto the line in the direction of V(f> is also a Cantor set of the same dimension for 
typical (full vo lum e _) M.M), see Q. T lmS the set of v alueS t „ lldmg 
to when L(u — t) is tangent to T form a Cantor set of dimension dim#(r). If tt is the 
projection from T onto a line in the direction of V</>, then the projected measure n^i>r has 
local dimension dimff(r) for typical (a, b, c, d). We have 



v{L (u — e)) = / / dm^udur. (36) 

J0 JL+(u-e) 

To estimate this integral we bound it above via the inequality m 1 u{^j u r\L + {u — e)) < C-y/e, 
and bound it below using the fact that for t > e/2, m 7 u(7 u n L + (u — e)) > Cy/e. Here 
C > is a uniform constant. Putting this together we obtain for typical (a, b, c, d) 

u(L + (u - e)) = f I dm^dur w yfe ■ e dimH ^ = e V2+dim H (r)_ (37) 

JO Jl+(u-e) 

D 



B. Numerical results 

We now examine the convergence of the numerically estimated values to the theoretically 
expected ones. For the numerical simulations, we rewrite observable ( 130]) in two forms 



06i,i, 06i,2 : M 3 — ;> M (form (130]) is recovered for suitable values of a, b, c, d): 

4> e ,i{x,y,z) = cos(2ir9)(x - x ) + sm(2ii9)(y - y ), (38) 

(f)e, 2 (x, y, z) = cos(27r6)(x - x ) + sin(27r#)(z - z Q ). (39) 

The level sets associated to ( 138]) and (139]) are planes orthogonal to (cos(27r#), sin(27r#), 0) 
and (cos(27r#), 0, sin(27r6>)), respectively. Figure [7] shows the dependence of the estimates 
of £ with respect to the block length N b i ock i en for both (138]) and ( 139]) at 9 = 0.5. Figure [7] 
suggests that the block length Nuockien = 10 4 is sufficient to get an estimate coherent with 
the theoretical value ( l34j) for observable (J38l) . whereas the same value is not sufficient for 
observable (l39|) . This is illustrated in Figure |8]for a range of values of 9: reliable estimation 
is obtained with block length Nuockien = 10 4 for observable ( 138]) (panel (A)) but not for 
observable ( 139]) (panel (B)), for which Nbiockien — 10 6 seems to suffice (panel (C)). 

In summary, the minimum block length required for (approximate) convergence to the 
theoretical value may vary strongly with the location within the attractor of the extremal 
point p in ( 132]) . that is with the relative position of the attractor and the level sets. Also, 
the minimum block length may depend on the dimensionality of the attractor. Numerical 
experiments suggest that reliable estimation is more difficult when the dimensionality of 
the attractor is smaller. Figure |9] indeed shows better agreement with the prediction 
of (134]) for the larger values of A which also correspond to a larger dimension according 
to (USD. 



Lastly, we consider observable ( 129]) . As we did in Figure [5] we illustrate the effect of the 
point pm "dropping out" of the attractor A. To achieve this, we iterate the solenoid map 
starting from an arbitrarily chosen initial condition. After discarding a transient of 10 5 
iterates, we regard the final point p° M of the orbit as being generic with respect to the 
Sinai-Ruelle-Bowen measure on A. Figure [10] (A) shows the sensitivity of the estimates 
of £ with respect to block length for observable ( 129]) when pu is equal to p° M as obtained 
above. The estimates display strong oscillations around the theoretical value and barely 
seem to settle for very large block lengths Nuockien > 10 7 . We return to this problem in 
Section |VJ 

We then choose pu as a perturbation of point p° M in the radial direction in IR 3 : namely we 
set pm = Pm = (1 + ^)Vm- By dissipativity of the solenoid map, we expect p l M $_ A with 
probability 1 when t ^ 0. We find out that when t is sufficiently large (Figure ITOlB)). the 
estimates of £ converge to the theoretically expected value (|34j) already for block lengths 
of 1000. However, when t is small (Figure [TU] C)) the estimates are closer to the value 



attained within the attractor (I3TJ) for small block lengths, whereas convergence to the 
theoretically expected value ( )34l) takes place for Nuockien larger than about 10 5 . 



V. NON-UNIFORMLY HYPERBOLIC EXAMPLES: THE HENON AND LOZI 
MAPS 

We here consider the Henon map |3j, |9( 

ha, b ■ K 2 ->• M 2 , h a , b {x, y) = (1 - ax 2 + y, bx), (40) 



for the classical parameter values (a, 6) = (1.4,0.3) and the Lozi map |12l |58| 

l a , b :R 2 ^R\ l ajh {x,y) = {l-a\x\+y,bx)i ( 41 ) 

for (a, b) = (1.7,0.1), under the observables 

<Pa(x, y) = - dist(p,p M )°, with p = (x, y) e M 2 , (42) 

<j>e(x,y) =xcos(2tt6) + ysm(2ir6), (43) 

where a > and 9 G [0, 2tt] are parameters and pm is a point in M. 2 . Following the 
discussion for the solenoid map, we could conjecture that 

1 dim^(A) 



f a 



for (f) = 4> a and £>m G A; (44) 



-- = dim^(A) - - for = <f> . (45) 

The numerical verification of these conjectures turns out to be rather problematic. First 
of all for a given system it may be very hard or even unfeasible to compute an estimate 
of the Hausdorff dimension. For this reason, we will use the Lyapunov (Kaplan- Yorke) 
dimension of the Henon or Lozi attractor instead of the Hausdorff dimension appearing 
in ( I45|) -( l4*4|) . The Lyapunov dimension of an attractor A C M n of a system with R n as 
phase space is defined as 

dim L (A) = k + ^i= lX \ (46) 

— Xfc+i 

where xi ^ X2 > • • • > Xn are the Lyapunov exponents and k is the maximum index for 
which Ylj=iXj ^ 0- It is believed that the Lyapunov dimension forms an upper bound 



for the Hausdorff dimension under general conditions 26|, |35 | 



For the Henon map under observable (TTBl) . and in view of the results of a recent pa- 
per J9(, it is expected that formula (j4"4"j) holds for so-called Benedicks- Carleson parameter 



values [23[. Such parameter values, however, are obtained by a perturbative argument 
near (a, b) = (2, 0), where the bound on the smallness of b is not explicit. Moreover, the 
parameter exclusion methods used to define the Benedicks- Carleson parameter values are 
not constructive. For these reasons, it is not possible to say whether Benedicks-Carleson 
behaviour is also attained at the "classical" parameter values (a, b) = (1.4,0.3). Despite 
this, (1441) forms our best guess for the value of £. 

For planar observables, we again study the tip set T C A as defined for the Solenoid map, 
namely, for fixed e > and p = (x, y), let 

r = {p e L + (u - e ) n A : T pl u (p) • V0(p) = 0}, (47) 

and for each e < eo, consider the set T(e) C T fl L + (u — e). Then for each p 6 T(e), there 
exists t < e such that 7 u (p) is tangent to L + (u — t). For the planar observable <fi we would 
expect to obtain (as with the solenoid): 

m f {7»nl + (u-e)} = 0(Vn), (48) 

where m 7 « is the conditional (Lebesgue) measure on the one- dimensional unstable mani- 
fold. However in this calculation we have assumed that the tangency between ^ u (p) and 
L{u — t) is quadratic, and that the unstable segment is sufficiently long so as to cross 
Liu — e) from end to end. For the Henon map both of these conditions can fail. In par- 
ticular, the Henon attractor admits a critical set of folds that correspond to points where 
the attractor curvature is large. More precisely the critical set is formed by homoclinic 
tangency points between stable and unstable manifolds. This set has zero measure, but 
it is dense in the attractor. Furthermore the attractor has complicated geometry, where 
local stable/unstable manifolds can fold back and forth upon themselves. However, the 
regions that correspond to these folds (of high curvature) occupy a set of small measure. 



See 57| and references therein for a more detailed discussion. 



To compute the tail index, we conjecture to have the following formula: 

- l - = dim^ (u) - l - (49) 

where v is the SRB measure for the Henon map (at Benedicks-Carleson parameters). This 
would follow from the estimate: 

u(L + (u -e))=[ f dm^udvr ~ v^ ■ e dim "( r \ (50) 

J0 JL+(u-e) 

where the factor of s/e comes from equation (|35|) . To obtain equation ( 1491 . we would 
need to show that dim^(r) = dim#(z/) — 1. This is perhaps harder to verify and it will 



depend on the regularity of the holonomy map taken along unstable leaves. Finally we 
would project this set onto a line in the direction of V0(jo), and typically the projection 
would preserve the dimension. 

Figure [UJ shows the dependence of the estimates of £ with respect to the block length 
Nbiockien for the Henon map under the observable ( Tl3l) . We see that the estimates exhibit 
strong oscillations around the value predicted by (144)) even for fairly large block lengths. 
Figure [12] shows the dependence of the estimates of £ with respect to the block length 
Nbiockien for the Henon map under the observables (fl3l) at 9 = 0, 0.5. The horizontal lines 
represent the values predicted by ( l4"5j) . where, as above, we have used the of Lyapunov 
(Kaplan- Yorke) dimension of the Henon attractor instead of the Hausdorff dimension. We 
see that block lengths of at least 10 4 are required for the estimation to reach the neigh- 
bourhood of the value predicted by ( ]4"5]) . However, the estimates still exhibit substantial 
oscillations around the predicted values for block lengths as large as 10 7 , although both 
the variability of the individual point estimates and the estimation uncertainty are here 
much less pronounced than in Figure [11] 

We had already seen the above behaviour in the solenoid map: namely, the estimates in 
panel (A) of Figure ITUl also exhibit larger variance and variability than those in panel (C). 
In that case, however, the theoretical value of panel (A) is not conjectural, since it follows 
from the theory discussed in Section IIVI for observables such as (129]) when the point pm 
belongs to the attractor. 

Hence we do not interpret the variability in Figures [JTJ and [12] as a dismissal of P5]) -( ]4^]) . 
Rather, we claim that this behaviour is due to a problematic aspect of the numerical 
estimation. To illustrate our claim, we more carefully examine the estimates of the GEV 
distribution obtained for block lengths of 5000 and 10000, for observable ( 1431) with 9 = 0. 
In this case, the observable simply coincides with the projection on the x-axis: this is very 
useful for the visualisation. 

The kernel-smoothed density of the block maxima show various peaks (panel (Al) in 
Figure IT3]) . A particularly pronounced peak occurs nearby x = 1.2727. Examination 
of the points on the time series of the block maxima (panel (Bl)) and of the points on 
the Henon attractor corresponding to the block maxima (panel (CI)) reveals that this 
peak is associated to a pair of branches of the attractor that exhibit a turning point 
slightly above 1.2727. This peak corresponds to a "corner" in the quantile-quantile plot 
(panel (Dl)) comparing the empirical distribution of the block maxima to fitted GEV 
distribution. For values of x at the left of the peak, the empirical distribution of the block 



maxima displays a strong deviation from the fitted GEV distribution. 

When the block length is increased to 10 4 (right column of Figure [T3]) . the kernel-smoothed 
density of the block maxima drops to almost zero at the left of the peak (panel (A2)). In- 
deed, the portion of the Henon attractor corresponding to the block maxima (panel (C2)) 
does no longer include the two leftmost branches which were found in panel (CI). More- 
over, a much smaller fraction of points now belongs to the branch of the attractor having 
a turning point at 1.2727. This also corresponds to the peak in the density being lower in 
panel (A2) than in panel (Al). More importantly, this correspond to a much better over- 
all fit to the GEV distribution: as illustrated by the quantile-quantile plot in panel (D2), 
there still is some deviation at the lower tail, but it is orders of magnitude smaller than 
in panel (Dl). 

We believe that this is the explanation for the poor convergence to the theoretical esti- 



mates which we have found in Figure [T2J also see 42j for a related discussion. The fractal 
structure portrayed in panels (Dl-2) of Figure [13] is indeed present at all spatial scales 
near the extremal point p = (x, y) on the Henon attractor for which observable (143p with 
9 = is maximised. As blocks of increasing lengths are used, increasingly many attrac- 
tor branches are discarded. Near the block length values for which one major branch is 
discarded, a better agreement is obtained between the sample of block maxima and the 
limiting GEV distribution. These are the block length values for which we expect the 
estimated value of £ to lie closer to the theoretical prediction in panel (A) of Figure [12] 

The effect of the variability in the estimates is illustrated in Figure HU where we show 
estimates of £ for observable ( |45|) with several values of 9 and with four block lengths. 
For Nuockien — 10 3 , the estimates vary substantially across the range of values of 9 (Fig- 
ure HU (A)). Varying 9 from to 1 amounts to rotate the level sets of the observable (145 p . 
which are straight lines. Hence, this amounts to slide the extremal point p for which 
observable (143]) is maximised on the Henon attractor (compare with (|7|)). The horizontal 
plateau in Figure [T4l (A), occurring approximately for 9 between [0.5,0.75], corresponds 
to the extremal point p belonging to the leftmost tip-like portions of the Henon attractor: 
large variations in 9 in this range correspond to small variations in p. 

For block lengths of N^i oc ki eri = 10 4 , (Figure [14] (B)), the estimates are more uniform 
across 9. The same holds for Nuockien = 10 5 and 10 6 and we see a definite bias in the 
latter case, which has the same sign and approximately the same value for all 9. Roughly 
speaking, choosing block lengths of at least Nuockien = 10 4 ensures that we only select 
block maxima in branches of the Henon attractor which are close to its outer "peel", 



compare with Figure [13] (C1-C2). However, this does not necessarily guarantee accurate 
estimation of the limit value of £, for the reason illustrated in Figure [TBI (D1-D2). 

We argue that the same explanation holds for the variability of the estimates in panel (A) 
of Figure [UJ and for the even poorer convergence in Figure [TTJ Plots similar to Figure [13] 
for the latter case suggest that as block length is increased, the probability mass that is 
lost at the lower tail of the empirical distribution of the block maxima is redistributed 
amongst other attractor branches which lie closer to the point p° M . To illustrate this 
process we chose observable (H3j) for ease of visualisation. 

Similar considerations hold for the Lozi map ()4ip . Figure [15] shows the sensitivity of 
the numerical estimates of £ with respect to the block length used to compute the 
maxima for observable fj42|) . For the chosen parameter values, we obtain the estimate 
dim^(A) = 1.185, in good agreement with the bounds 1.176669 < dim#(A) < 1.247848 
on the Hausdorff dimension of the Lozi attractor A proved in [34|. When the point pm is 
chosen in the attractor of the Lozi map, the estimates display strong oscillations around 
the value predicted by the theory (Figure [151 panel (A)), as in Figure ITO We then choose 
Pm — (0.2, 0.01): this point does not lie on the attractor of the Lozi map, but the nearest 
point (x, y) on the attractor belongs to one of the straight portions. Also in this case we 
observe oscillations around the theoretically expected value (Figure [TBI panel (B)), like 
in Figure | 



VI. THE LORENZ63 AND LORENZ84 FLOWS 



The theoretical and numerical machinery developed in the previous sections is now ap- 



plied to two paradigmatic ordinary differentia 
Lorenz. We first of all consider the model of 



equations, both derived and studied by Ed 



38 



x = a(y — x),r 

y = x(p-z)-y, (51) 

z = xy — (3z, 

derived from the Rayleigh equations for convection in a fluid layer between two plates. 
Here a is the Prandtl and p the Rayleigh number. We refer to this as the Lorenz63 model 
and fix a = 10, = 8/3 and p = 28, which is a fairly common choice in the vast literature 



,fl 



on the Lorenz63 system, see e.g. |l|, |46|, |48j . The statistics of extremes has been previously 



analysed in 52( , who found smooth-like variation of the GEV parameters with respect to 



changes in the parameter p within a suitable range. 



We also study a three-dimensional system proposed by Lorenz in 1984 39| : 



X = 


—ax — y 2 — z 2 + aF. 


y = 


—y + xy — bxz + G, 


z = 


—z + bxy + xz. 



(52) 



This is derived by a Galerkin projection from an infinite dimensional model for the at- 
mospheric circulation at mid-latitudes in the Northern Hemisphere. The variable x is the 
strength of the symmetric, globally averaged westerly wind current. The variables y and 
z are the strength of cosine and sine phases of a chain of superposed waves transporting 
heat poleward. The terms in h represent displacement of the waves due to interaction with 
the westerly wind. The coefficient a, if less than one, allows the westerly wind current 
to damp less rapidly than the waves. The time scale of t corresponds to about 5 days. 
The terms in F and G are thermal forcings: F represents the symmetric cross-latitude 
heating contrast and G accounts for the asymmetric heating contrast between oceans and 
continents. System fl52l) has been used in both climatological studies, for example by 



Sever a 
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M 



works have 



49j. Almost 



coupling it with a low- dimensional model for ocean dynamics [50 [ 

examined its bifurcations, mainly in the (F, G)-parameter plane 

nothing is known theoretically regarding the structure of its strange attractors. As in the 

above references, we fix a = 0.25 and b = 4 and consider the chaotic dynamics occurring 

at (F, GO = (8,1). 

We analyse time series generated by observables computed along orbits of these flows, 
sampled every At time units. We fix At = 0.05 time units for the Lorenz63 and At = 0.1 
for the Lorenz84 model. We consider the two observables 

<fri(x,y) = - dist (p,p M ), with p = (x, y, z) € M 3 , (53) 

</> 2 (x,y,z)=x. (54) 



Observable 02 has a clear physical meaning in both models: for ( jSTjl . the variable x repre- 
sents the intensity of the convection, whereas in ( )52l) the variable x represents the strength 
of the westerly wind current. As in the previous sections, we examine the sensitivity of the 
numerical estimates of £ with respect to the block length used to compute the maxima. 

We first consider the Lorenz63 system (15TT) . It will be useful to recall some geometrical 
facts of the Poincare map to z — constant sections. Given the planar sections S = 
{(x, y, 1) : |x|, \y\ < 1}, and £' = {(1, y, z) : \y\, \z\ < 1}, the map P : S — )■ E decomposes 



as P = Pi o P l7 where P x : E — > E' and P 2 : E' — > E. To describe the form of P, let 
/3 = |A S |/A U , /?' = |A SS |/A U , where A s , \ ss and A„ are the eigenvalues of the linearised 
Lorenz63 flow at the origin, with X s = —8/3, X ss = —22.83 and A M = 11.83 for our 
choice of parameters. Then it can be shown that P\(x,y, 1) = (l,x@ y,x@), and P2 is a 



diffeomorphism, see [32j. Thus the rectangle E + = {(x, y, 1) : x > 0, \y\ < 1} gets mapped 
into a region Pi(E + ) with a cusp at y — 0. The cusp boundary can be represented as the 
graph \y\ = z^ '& ~ z 8 . The flow has a strong stable foliation, and we form the quotient 
space E = E/ ~ by defining an equivalence relation p ~ q if p e / ~f s {q), for a stable leave 
7 s . Hence the map P : E — > E can be reduced to a uniformly expanding one-dimensional 
map / : E — > E, with a derivative singularity at a; = 0. Here E identified with [—1, 1], 
and f'(x) ~ jxp -1 near x = 0. 

The Lorenz flow admits an SRB measure z/ which can be written as v = vp x Leb (up to a 
normalisation constant). The measure vp is the SRB measure associated to the Poincare 
map P, and the measure is exact dimensional, i.e. the local dimension is defined z/-a.e., 
see 25j . Using the existence of the stable foliation, and the SRB property of u, we can 
write up as the (local) product u^u x u^s where u^ is the conditional measure on unstable 
manifolds, and u^ is the conditional measure on stable manifolds. We can identify each 
measure i/ 7 « (via a holonomy map) with that of the invariant measure Uf associated to 
/. The measure uj is absolutely continuous with respect to Lebesgue measure, but it has 
zero density at the endpoints of E, that is 

u f ([l - e, 1]) w e 1//3 w e 4 ' 4 as e -»■ 0. (55) 

From this analysis we can now conjecture the values of £. Following the reasoning as 
applied in Section [IV] the conjectural values of £ are 

— - = dimji(u), for observable ( J53l) . (56) 

111- 

— = — H hG?s where rf s <l for observable (154|) . (57) 

The constant <i s comes from the dimension of u s which is (numerically) seen to be small 
due to the strong stable foliation. As in Section [V] we replace the Hausdorff dimension 
with the Lyapunov dimension, which we numerically estimate at dim^A m 2.06. We 
take this value to be the estimate of the local dimension of u. In contrast with the 
solenoid and Henon maps, the tail index associated to observable ( 153]) comes from an 
estimate of the measure of u(L + (u — e)) which we assume scales as the product of the 
three factors: y/e ■ e du ■ e 8da . Here the factor \fe comes from the measure u conditioned 



on A n L + (u — e) in the (central)-flow direction, while the factor e du comes from the up- 
measure conditioned on unstable manifolds that terminate at the cusp. In a generic case 
we would expect d u — 1. However, since we are near the cusp (namely near the boundary 
<9X) we have d u = 1/(3 = 4.4 due to the zero in the density of Uf, see ( )55|) . Finally we 
have a contributing factor e 8da that comes from the the strength of the cusp at P(<9£), 
with d s the local dimension of i/ 7 ». We would expect typically that d s ~ 0.06, but it could 
be much smaller if we are in the vicinity of the cusp. 

For the numerical simulations we first consider observable (|53l) . where pm is a point chosen 
in the attractor by the same procedure used before, namely selecting the final point of 
an orbit of length 10 3 time units. The estimates converge to the theoretically expected 
values of — 1/dimL A rs —0.5, see Figure [TBI (A). Note that convergence is attained here 
for block lengths of a few thousands, unlikely what has been observed for the Henon and 
Lozi maps. We also obtain convergence to the conjectured value (|5"T|) for observable (154|) . 
see Figure [TBI (B). Also in this case the convergence is much faster than for the Henon 
and Lozi maps. 

For the Lorenz84 system ([521) . Lorenz detected a Henon like structure in a Poincare 



section with the plane y — 0, see [39j, Figures 7 and 8]. If this conjectural structure was 
correct, then the attractor would coincide with the two-dimensional unstable manifold of 
a saddle-like periodic orbit of the flow of ( ]52]) . 



Assuming that there is exists an SRB measure v supported on this attractor, and that 
there is a local product structure so that v can be written as v^u x v^» (as with Lorenz63), 
then following the reasoning of Section ITVl the conjectural values of £ are 

— - = dim#(i>), for observable ([53]) . (58) 



1 dmif^^y 



+ diniff(i/ 7 s) for observable (154|) . with dim(7 u ) = 2. (59) 



£ 2 

In this conjecture, it is assumed that the level sets L + (u — e) meet the unstable manifolds 
via generic quadratic tangencies (unlike Lorenz63). The estimates for observable ( l53l) 
display oscillations around the theoretical value (|58|) . see Figure [T7] (A). This behaviour 
similar to what observed for the Henon and Lozi maps, see Figs. [11] and [15] (A). The 
estimates in Figure IT7l(B) display oscillations around the value ( 159]) : again this is similar 
to what was observed in the Henon and Lozi maps, see Figs. [12] and [15] (B). 



VII. DISCUSSION AND CONCLUSIONS 

This paper has presented an extension of the currently available extreme value theory 
for dynamical systems to types of observables which are more similar to those found 
in applications. Namely, the observables considered here are not (necessarily) functions 
of the distance from a point which is generic with respect to the invariant measure of 
the chaotic system. Formula ( 134)) and its generalisation fl59|) were derived under generic 
assumptions on the geometry of the invariant manifolds underlying the strange attractor. 
Current research by the Authors aims at formulating explicit conditions under which such 
formulas hold, both for uniformly and non-uniformly hyperbolic systems. Preliminary 
findings suggest the following. Suppose we have a system with an attractor A C M. d 
that supports a Sinai-Ruelle-Bowen (SRB) measure v. Moreover suppose that A admits a 
local product structure so that v can be locally regarded as the product measure u 7 u x Vy, 
where v^ (resp. v^s) are the conditional measures on unstable (resp. stable) manifolds. 
Since v is SRB, the measures v^v. are equivalent to the Riemannian measures on the 
unstable manifolds, and their local dimension d u is an integer. The local dimension of 
u-ys is typically non-integer. For sufficiently smooth observables (p : M. d — > K. that have 
maxima off A, we conjecture that the tail index £ is given by the formula: 

The factor 4^ comes from assuming that the level sets meet the unstable manifolds in 
generic (quadratic) tangencies. The factor d s is the local dimension of Vy. We believe 
that this dimension d s is equal to the dimension of the tip set V as defined by equation 
( |33l) . Most of our examples had dim#(r) < 1, but in general this could be larger than 
1. If this is so, then the projection of V onto a line in the direction of V0(p) would 
typically have dimension equal to one, and the intersection of V with each level set would 
(typically) be an uncountable set of positive Hausdorff dimension. Thus in addition to 
studying regularity of unstable holonomies, a careful analysis of the attractor's geometry 
would be required when estimating the //-measure of the level regions nearby the extremal 
point p. 



It is of interest to verify the above formula for maps where d u is 
is the case for the so-called quasi-periodic Henon-like attractors 



arger than one: such 
54 l55| . which are 



contained in the closure of the 2D unstable manifold of a saddle-like invariant circle. For 



flows, this situation corresponds to d u = 3, see e.g. 



4J]. Also, the Lorenz63 example 



presented in Section [VI] shows beyond doubt that the geometry of the attractor can play 



a substantial role in determining the limit GEV distribution. In that case the level sets of 
the observable do not meet the attractor via quadratic tangencies: instead, the level sets 
meet the attractor at cusps where the measure v^u has a zero. Therefore relation ( 1601 
fails to hold and the alternative formula (1571) is derived. This situation bears resemblance 
to the configuration Figure [11(C) for Thorn's map, which leads to formula (ITT)) for the tail 
index. Similarly, a modified formula for £ is expected to hold for the Lozi map under the 
observable 4>(x, y) = x, for which the extremal point p coincides with a cusp-like point in 
the attractor. 

As far as applications are concerned, this paper both points at the further development 
of useful methodologies and also raises a number of significant questions. We envisage 
the development of estimation methods for the parameters of the GEV distribution which 
take into account the information provided by formulas such as (160)) . Given a concrete 
system, parameter estimation would be complemented by an analysis of the structure of 
the attractor to determine appropriate values for d u and d s . Specifically, d u could be 
estimated by examining Poincare sections of the attractor and/or finite time Lyapunov 
exponents. Calculation of Lyapunov exponents would then yield d s through the relation 
dime (A) = d u +d s , which follows from the local product structure of the invariant measure. 
Such analysis would also aim to ascertain whether a formula like ( 160)) or appropriate 
modifications like ( )5Tj) should be used. This information could be fed into the parameter 
estimation procedure in an appropriate Bayesian setting. 

In the presence of parameter-dependent systems, these formulas provide an explanation 
for the smooth-like dependence of extreme value statistics with respect to changes in the 



control parameters. This phenomenon was first observed in la, [l7| and the implications 



for parameter estimation in non- stationary systems were discussed in 52j. This phe- 
nomenon critically depends on the structure of the observables: indeed, for observables 
like (|5]) we expect formulas like (160)) or ( )5Tj) . which could depend smoothly on control 
parameters through smooth-like dependence of the attractor dimension on control param- 
eters. On the other hand such smooth-like dependence is rather unlikely to occur for the 
observables considered so far in the theoretical work, which are of the form (j4)). Indeed, 
SRB measures in geophysical systems are usually singular with the Lebesgue measure in 
phase space, due to dissipation. Therefore, even if the point p M is generic for the SRB 
measure for a given value of the control parameters, this situation is typically not stable 
under parameter variation. If pu is fixed, then one would expect jumps in the value of 
£ whenever pu "drops off" or "drops into" the attractor, see the discussion for Figures [6] 



and US 

As far as the questions are concerned, the main one appears to be the extremely 
slow convergence displayed by the Henon-like attractors considered here (see e.g. Fig- 
ures [HI [T5l[T7l) . Such a slow convergence has been previously observed in more complex 
atmospheric models, see [51 |. Does such a slow convergence take place in state-of-the-art 
global climate models? This might pose a very serious methodological problem for those 
studies aiming at quantifying climatic change in extremes, for example changes in the 
behaviour of hurricanes, wind storms and extreme rainfall. 

These problems even raise the following provocative question: how relevant are limit 
laws for extreme behaviour, if it takes too long for the limit to be attained for any 
practical purpose? This question may have different answers. One possibility is that 
novel modelling approaches could be developed to provide more reliable estimates of 
extreme behaviour, not necessarily restricted to the standard limit laws such as the GEV 



or the Generalised Pareto distributions lOf. Alternatively, novel parameter estimation 



procedures might be developed, that incorporate corrections or modifications to account 
for the phenomena illustrated for Henon-like attractors, also see Figure [13j The results of 
this paper seem to suggest that whatever the final answer (s), the methods will have to take 
into account the geometry and the fractal nature of the strange attractors underlying the 
dynamics. We believe that these questions and problems will be the subject of significant 
research efforts in the near future. 



Appendix A: Parameter estimation for the GEV distribution 

We now describe the procedure which we have used to estimate the parameters /i, a, £ of 
the GEV distribution (J2j). Consider Nb max values x\, . . . , %N hmax which we assume to form 
a random sample from ([2]). For the systems under consideration, it often turns out that 
the theoretically expected value of £ is smaller than —0.5. In such cases, the standard 
maximum likelihood approach cannot be used, because the maximum likelihood estimator 



is not regular 10j. We therefore resort to the method of L-moments 



331. For the GEV 



distribution, the L-moments estimation equations are 



see Table 1 in 



Ai = /x-|(l-r(l-0), (Al) 

A 2 = -|(l-2«)r(l-0, (A2) 



33| . Given the sample x±, . . . ,%N bmax , we use the R package Lmoments 
(http://cran.r-project.org/) to estimate the first three L-moments Aj, i = 1,2,3. 
Eq. (1A3h is then solved for £ by a Newton method, starting from the initial estimate 



£ = 7.859^ + 2.9554z 2 , with 2 = 2/(3 + jj) - log 2/ log 3, see Table 2 in |33|. Once an 
estimate of £ is obtained, this is plugged into ( 1A2I) . which is solved for er. Lastly (lAljl is 
solved for /x. 

For the numerical computations, which also include quantifying the estimation uncer- 
tainty, we adopt the following procedure. Positive integers Nb max , Nuockien and N samp are 
fixed: here Nbmax is the total number of block maxima to be computed and Nbi 0C Hen is 
the length of the data blocks over which each maximum has to be extracted. We first 
discard a transient of 10 5 iterates with the map under consideration (e.g. the Henon or the 
solenoid map). Then a total of Nbmax • Nuockien iterates with the map is computed and the 
Nbmax block maxima are extracted. This sample is divided into N samp sub-samples, each 
containing Nbmax /N samp values. We then apply the above L-moment estimation procedure 
to each sub-sample, thereby obtaining N samp distinct parameter estimates 

{(/i s ,a s ,£ s ) \s = l,...,N samp }. (A4) 

The sample means of the estimates ( 1A4I) : 

N aamv W.„™„ AT„„_ 



/' 



J2 ^, = — — Yl °» i = m — E & ( A5 ) 



1 ' samp -, ly samp -. ly samp -, 



are taken as the final GEV parameter estimates and the standard deviations 

^ iV samp , -IN samp , IN samp 

4 = jt- E ^-^ 4 = y- E (—-) 2 ' 4 = t- E fe-o 2 ( A6 ) 

iv samp -, iv samp -, ly samp -, 



are taken as estimates of uncertainty for the final values ( 1A5I) . 
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Figure 1: Sketch of the three situations considered in Theorem IIII.ll for the level sets L(u) 
(defined in (J7J)) for observable cp a (fT3|) . 
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Figure 2: Sketch of a few possible configurations for the level sets L(u) (defined in ([7])) for 
observable (j) ab flU}. (A) (a, 6) = (2,1.25); (B) (a, 6) = (0.5,0.75); (C) (a, b) = (1.5,0.7). 
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Figure 3: Point estimates (crosses) and estimation uncertainty (vertical bars) of the tail index 
£ versus block length N^i^k^n for Thorn's map (|12p under the observable (|14p with a = 2 
and p M = (0.510001,0.5090001) fixed, where 6/2 = 0.5 (left) and 6/2 = 1.75 (right). The 
horizontal dashed lines represent theoretically expected values according to (|18p . Crosses and 



vertical bars are the mean and ± one standard deviation of a sample of N, 



samp 



100 individual 



estimates along a single orbit. Individual estimates are obtained by the method of L-moments 
with sequences of Nb max = 50000 block maxima over blocks of length Nu oc ki en , as described in 
Appendix El 
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Figure 4: Point estimates (crosses) and estimation uncertainty (vertical bars) of the tail in- 
dex £ versus parameter b for Thorn's map (|12h under the observable (I14j) with a = 2 and 
p M = (0.510001,0.5090001) fixed and varying 6/2 = 0.01,0.1,0.25,0.5,0.75,1,1.25,1.5,1.75,2. 
The dashed line represents theoretically expected values according to (|18p . Point and interval 
estimates are obtained by the method of L-moments as for Figure HI with N\ miax = 10000, 



Nuockien = 10000 and N, 



samp 



100, see Appendix 1X1 
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Figure 5: Point estimates (crosses) and estimation uncertainty (vertical bars) of the tail index 
£ versus block length N\,i oc ki en for Thorn's map Q12|) under the observable (|13p with a = 1 and 
Pm = (xm,Vm) with 2/m = 0.510001 fixed and (A) xm = 0.9, (B) xm = 10, (C) xm = 1-01, 
(D) xm = 1-1. Horizontal lines labelled by 1, 2 represent theoretical values according to Q15|) 
and (|16h , respectively. Estimates are obtained by the method of L-moments as for Figure HJ 
with Nb max = 50000 and N samp = 100, see Appendix lAl 
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Figure 6: Point estimates (crosses) and estimation uncertainty (vertical bars) of the tail index 
£ for Thorn's map (|12|) under the observable (|13|) with a = 2 and (A) yu = 0.5090001 and 
xm varying, and (B) xm = 1.2 and yu varying. Horizontal lines labelled by 1, 2, 3 represent 
theoretical values according to (|15p , (|16|) and (|17|) , respectively. The method of L-moments was 



used as described in Appendix [7\] with N bmax = 10000 and N. 
10 4 (red) and N blocklen = 10 5 (blue). 
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Figure 7: Point estimates (crosses) and estimation uncertainty (vertical bars) of the tail in- 
dex £ versus block length Nbi oc ki en for the Solenoid map (f2"T|) under the observable ([38]) (left) 
and (|39|) (right) with 9 = 0.5. The horizontal dashed lines represent theoretically expected 
values according to (|34p Estimates are obtained by the method of L-moments as for Figure [U 



with Nto 



10000 and N, 
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100, see Appendix lAl 
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Figure 8: Point estimates (crosses) and estimation uncertainty (vertical bars) of the tail index 
£ versus parameter 9 for the solenoid map (j27|) under the observable (j38|) (A) and (f39|) (B,C) 



with xq = yQ = zq = 3, where 9 



10 



for i = 0, . . . , 9. The dashed line represents theoretically 



expected values according to (|34p . Estimates are obtained by the method of L-moments as for 



Figure [U with N^ 



10000 and N. 
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blocklen 



10 6 (B). 
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100, see Appendix|A3 where Nm 
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Figure 9: Point estimates (crosses) and estimation uncertainty (vertical bars) of the tail index 
£ versus parameter A for the solenoid map (I27|) under the observable (j38|) with 9 = xq = yo = 
0. The dashed line represents theoretically expected values according to (|34p . Estimates are 
obtained by the method of L-moments as for Figure [U with N}y max = 50000, N samp = 100 and 
Nbiockien = 10000, see Appendix [XJ 
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Figure 10: Point estimates (crosses) and estimation uncertainty (vertical bars) of the tail index 
£ versus block length Nuocklen for the solenoid map (|27h under the observable (|29h with a = 0.3. 
(A) pm is chosen as a point p° M G A as described in the text. (B) j»m = V 1 m = (1 + *)Pm with 
t = 0.1. (C) pjvf = p f M with £ = 1. The dashed lines represent theoretically expected values 
according to (I15p Estimates are obtained by the method of L-moments as for Figure HI with 



Nbmax = 10000 and JV, 



samp 



100, see Appendix lAl 
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Figure 11: Point estimates (crosses) and estimation uncertainty (vertical bars) of the tail index 
£ versus block length A^/ oc fcz en for the Henon map (|10"j) under the observable (pi2"j) with a = 2. 
The horizontal dashed line represents theoretically expected values according to (|44p . with the 
Lyapunov dimension replacing the Hausdorff dimension, see text. Estimates are obtained by the 



method of L-moments as for Figure HI with Nj, n 



50000 and N 



samp 



100, see Appendix lAl 
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Figure 12: Point estimates (crosses) and estimation uncertainty (vertical bars) of the tail index 
£ versus block length A^ oc fc/ en for the Henon map (pTOj) under the observable (|35]l with # = 
(left) and 9 = 0.5 (right). The horizontal dashed lines represent theoretically expected values 
according to (|45[) . Estimates are obtained by the method of L-moments as for Figure HI with 



N, 



brnax 



10000, and N, 



samp 



100, see Appendix |A"1 
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Figure 13: Diagnostics of the GEV distribution fit for the Henon map (|40p under observable (|43h 
with 9 = 0, with Ni> max = 5 • 10 5 block maxima computed over blocks of length Nbiockien — 5 • 10 
(left column, Al-Dl) and Nbiockien — 1-2- 10 5 (right column A2-D2). (Al,2) Non-parametric log- 
densities of the block maxima, obtained by Gaussian kernel smoothing with bandwidth 0.000002. 
(Bl,2) Time series of the 5 • 10 block maxima (with the block sequential index on the vertical 
axis). (CI, 2) Points on the Henon attractor corresponding to the block maxima used in (Al) 
and (A2), respectively. (Dl,2) Quantile-quantile plot of the empirical distribution of the block 
maxima (horizontally) versus the fitted GEV distribution. 



U/ 



-0.6 
-0.7 
-0.8 
-0.9 
-I 
-1.1 
-1.2 
-1.3 
-1.4 



estimated 
theoretical 



(A) 



-0.6 rr 
-0.7 - 
-0.8 
-0.9 - 

-1 
-1.1 
-1.2 
-1.3 : 
-1.4 L 



estimated 
theoretical 



(B) 



0.25 



0.5 



0.75 







0.25 



0.5 



0.75 



-0.6 
-0.7 
-0.8 
-0.9 
-1 
-1.1 
-1.2 
-1.3 
-1.4 



estimated 
theoretical 



(Q 



-0.6 rr 
-0.7 - 
-0.8 
-0.9 - 

-1 
-1.1 
-1.2 P- 
-1.3 
-1.4 



estimated 
theoretical 



(D) 



0.25 



0.5 



0.75 







0.25 



0.5 



0.75 



Figure 14: Point estimates (crosses) and estimation uncertainty (vertical bars) of the tail index 
£ versus 6 for the Henon map (j40j) under the observable (J4~3l) for block lengths of (A) Nuockien = 
10 3 , (B) Nuocklen = !0 4 , (C) Nuoeklen = 10 5 , (D) N block i en = 10 6 . The horizontal dashed lines 
represent theoretically expected values according to (|45p . Estimates are obtained by the method 



of L-moments as for Figure IU with Nt> Ti 



10 5 , and N. 



samp 



100, see Appendix lAl 
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Figure 15: Point estimates (crosses) and estimation uncertainty (vertical bars) of the tail index 
£ versus block length Nbi oc ki e n f° r the Lozi map (|4U|) under the observable (|42|) with a = 2 and 
(A) for a point pm belonging to the attractor; (B) for pu = (0.2,0.01) The horizontal dashed 
line represents theoretically expected values according to (|44p . with the Lyapunov dimension 
replacing the Hausdorff dimension, see text. Estimates are obtained by the method of L-moments 



as for Figure m with N bmax = 50000 and iV< 



amp 



100, see Appendix lAl 
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Figure 16: Point estimates (crosses) and estimation uncertainty (vertical bars) of the tail index 
£ versus block length Nbi oc ki e n for he Lorenz63 flow (f5T|) (A) under the observable ([53]) where 
Pm is chosen as the final point of an orbit of length 10 3 time units starting from an arbitrary 
point, and (B) under observable (I54p . Horizontal lines labelled by represent theoretical values. 
Estimates are obtained by the method of L-moments as for Figure HI with Nf, max = 20000 and 



Ngamp = 100, see Appendix [Al 



-0.1 
-0.2 - 
-0.3 - 
-0.4 - 
-0.5 
-0.6 
-0.7 
-0.8 



(A) 

100 



1000 



L-mom estimate 

L-mom estimate 

theoretical 



10000 100000 

-'" blocklen 



**w 



le+06 



-0.5 | - 

-0.6 

-0.7 - - 

-0.8 

-0.9 

-1 - 

-1.1 

-1.2 

-1.3 -(B) 

-1.4 - 

100 



1000 



L-mom estimate ^ 
L-mom estimate ^ 
theoretical 

10000 100000 

1* blocklen 



le+06 



Figure 17: Point estimates (crosses) and estimation uncertainty (vertical bars) of the tail index 
£ versus block length Nuocklen f° r ne Lorenz63 flow (|5Tj) (A) under the observable ([55]) where 
Pm is chosen as the final point of an orbit of length 10 3 time units starting from an arbitrary 
point, and (B) under observable (|54p . Horizontal lines labelled by represent theoretical values 
according to (|58p (A) and (|59p (B). Estimates are obtained by the method of L-moments as for 
Figured with N bmax = 20000 and N samp = 100, see Appendix [A] 



